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Using the Trace Minimization Algorithm, we carried out an exact calculation of entanglement in 
a 19-site two-dimensional transverse Ising model. This model consists of a set of localized spin-i 
particles in a two dimensional triangular lattice coupled through exchange interaction J and subject 
to an external magnetic field of strength h. We demonstrate, for such a class of two-dimensional 
magnetic systems, that entanglement can be controlled and tuned by varying the parameter A = h/ J 
in the Hamiltonian and by introducing impurities into the systems. Examining the derivative of the 
concurrence as a function of A shows that the system exhibits a quantum phase transition at about 

(N : 

Ac = 3.01, a transition induced by quantum fluctuations at the absolute zero of temperature. 
I PACS numbers; 03.67.Mn 



I. INTRODUCTION 



, Entanglement, which is a quantum mechanical property that has no classical analog, has been viewed as a resource 

O I of quantum information and computation [H-Q- Intensive researches of entanglement measurement, entanglement 

^ monotone, criteria for distinguishing separable from entangled pure states and all the extensions from bipartite to 

^ • multipartite systems have been carried out both qualitatively and quantitatively [l| . At the interface between quantum 

' information and statistical mechanics, there has been particular analysis of entanglement in quantum critical models 

The critical properties in the entanglement allow for a screening of the qualitative change of the system and a deeper 
characterization of the ground state wavefunction undergoing a phase transition. At T=0, ground states of many-body 
■ systems contain all correlations concerning phases of matters. Traditionally, systems have been studied by looking, for 
^ example, at their external perturbations, various order parameters and excitation spectrum Methods developed 
from quantum information shed light on new ways of studying many-body systems |13l4l6l | , such as providing support 
for numerical calculations, like density matrix renormalization or design of new efficient simulation strategies for 
many-body systems. 

' Entanglement close to quantum phase transitions was originally analyzed by Osborne and Nielsen [l^ , and Osterloh 
. et al. |9|] for the Ising model in one dimension. Recently, we studied a set of localized spins coupled through 

Qv^ ' exchange interaction and subject to an external magnetic filed [l7l - l20| . We demonstrated for such a class of one- 
, dimensional magnetic systems, that entanglement can be controlled and tuned by varying the anisotropy parameter 
IL* ' in the Hamiltonian and by introducing impurities into the systems. In particular, under certain conditions, the 

.5^ ' entanglement is zero up to a critical point Ac, where a quantum phase transition occurs, and is different from zero 

X' 



above Ac [21 1 



?H In two and higher dimensions nearly all calculations for spin systems were obtained by means of numerical simula- 

9^ . tions [2^ [2^ . The concurrence and localizable entanglement in two-dimensional quantum XY and XXZ models were 
considered using quantum Monte Carlo [U, [2^. The results of these calculations were qualitatively similar to the 
one-dimensional case, but entanglement is much smaller in magnitude. Moreover, the maximum in the concurrence 
occurs at a position closer to the critical point than in the one-dimensional case [l| . 

The Trace Minimization Algorithm for Hermitian eigenvalue problems, like those under consideration in this paper, 
was introduced in 1982 by A. Samch and J. Wisniewski [2^. It was designed specifically to handle very large problems 
on parallel computing platforms for obtaining the smallest eigenpairs. Later, a similar algorithm (Jacobi-Davidson) 
for the same eigenvalue problem was introduced by Sleijpen and Van der Vorst in 1996. A comparison of the two 
schemes by A. Sameh and Z. Tong in 2000 [131 showed that the Trace Minimization scheme is more robust and efficient 
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In this paper, we use the Trace Minimization Algorithm [26|,|27| to carry out an exact calculation of entanglement in 
a 19-site two dimensional transverse Ising model. We classify the ground state properties according to its entanglement 
for certain class on two-dimensional magnetic systems and demonstrate that entanglement can be controlled and tuned 
by varying the parameter A = h/ J in the Hamiltonian and by introducing impurities into the systems. The paper is 
organized as follows. In the next section, we give a brief overview of the model, entanglement of formation and the 
trace minimization algorithm. Detailed methods are addressed in the appendix. We then proceed with the results and 
discussions of 1. the calculation of exact entanglement of a 19-site spin system, 2. the relationship of entanglement 
and quantum phase transition, and 3. the effects of impurities on the entanglement. The conclusions and the outlook 
are presented in the concluding section. 



II. METHOD 



A. Model 



We consider a set of localized spin-i particles in a two dimensional triangular lattice coupled through exchange 
interaction J and subject to an external magnetic field of strength h. The Hamiltonian for such a system is given by 

H = -Y,J,^,afa-^~hY,at, (1) 

<i,j> i 

where < i, j > is a pair of nearest-neighbors sites on the lattice, Jij = J for all sites except the sites nearest to the 
impurity site im, while around the impurity site Jij = (1 -f a)J, a measures the strength of the impurity which is 
located at site im, and af and erf are the Pauli matrices. For this model it is convenient to define a dimensionless 
coupling constant A = h/J. 



B. Entanglement of formation 

We confine our interest to the entanglement of two spins, at any position i and j |9j]. We adopt the entanglement of 
formation, a well known measure of entanglement [1^, to quantify our entanglement 21 1. All the information needed 
in this case is contained in the reduced density matrix pij. Wootters [29j has shown, for a pair of binary qubits, that 
the concurrence C, which goes from to 1, can be taken as a measure of entanglement. The concurrence between 
sites i and j is defined as [2^ 

C(p) = max{0, ei - 62 - €3 - £4}, (2) 

where the e^'s are the eigenvalues of the Hermitian matrix R = ^J~^fp(Q^ with p — {a^ ® cr^)p*(cr^ ® cr^) and is 
the Pauli matrix of the spin in y direction. For a pair of qubits the entanglement can be written as, 

E{p) = (3) 

where e is a function of the "concurrence" C 



where h is the binary entropy function 

h{x) = —x\og2X — {\ — x)log2(\ — x). (5) 

In this case, the entanglement of formation is given in terms of another entanglement measure, the concurrence C. 
The entanglement of formation varies monotonically with the concurrence. 

The matrix elements of the reduced density matrix needed for calculating the concurrence are obtained numerically 
using the formalism developed in the following section. 
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C. Trace minimization algorithm (Tracemin) 

Diagonalizing a 2^^ by 2^^ Hamiltonian matrix and partially tracing its density matrix is a numerically difficult 
task. We propose to compute the entanglement of formation, first by applying the trace minimization algorithm 
(Tracemin) (26l . [27j to obtain the eigenvalues and eigenvectors of the constructed Hamiltonian. Then, we use these 
eigenpairs and new techniques detailed in the appendix to build partially traced density matrix. 

The trace minimization algorithm was developed for computing a few of the smallest eigenvalues and the corre- 
sponding eigenvectors of the large sparse generalized eigenvalue problem 

AU = BUY., (6) 

where matrices A, B G ^nx" are Hermitian positive definite, U = [ui,...,Up] € C"^^ and S G MP^p is a diagonal 
matrix. The main idea of Tracemin is that minimizing Tr{X* AX), subject to the constraints X*BX = I, is equivalent 
to finding the eigenvectors U corresponding to the p smallest eigenvalues. This consequence of Courant-Fischer 
Theorem can be restated as 

^rmn^Tr{X*AX) = Tr{U*AU) = A„ (7) 

1=1 

where / is the identity matrix. The following steps constitute a single iteration of the Tracemin algorithm: 

• G = XlBXk (compute G) 

• G = VilV* (compute the spectral decomposition of G) 

• H = Q*AQ (compute H, where Q = Xfe^-i/^) 

• H ^ Wil.W* (compute the spectral decomposition of H) 

• Xk^ QW (now X^AXk = A and X^BX^ = I) 

• Xk+i = Xk-D ( D is determined s.t. Tr{Xl^^AXk+i) < Tr{X*AXk) ). 

In order to find the optimal update D in the last step, we enforce the natural constraint X^BD = 0, and obtain 



A BXk\ f D\ AX. 



XIB I \ Ll \ 



(8) 



Considering the orthogonal projector P ~ BXkiX^B'^Xk) ^X^B and letting D ~ {I - P)D, the hnear system ^ 
can be rewritten in the following form 

[I - P)A{I - P)D = [I - P)AXk. (9) 

Notice that the Conjugate Gradient method can be used to solve ([5]), since it can be shown that the residual and 
search directions r,p g Range(P)^. Also, notice that the linear system ([9]) need to be solved only to a fixed relative 
precision at every iteration of Tracemin. 

A reduced density matrix, built from the ground state which is obtained by Tracemin, is usually constructed as 
follows: diagonalize the system Hamiltonian H{X), retrieve the ground state j^" > as a function of A = h/J, build 
the density matrix p ~ j^" >< ^'|, and trace out contributions of all the other spins in density matrix to get reduced 
density matrix by p{i,j) = J2p < Ui{A)\ < Vp{B)\p\uj{A) > \vp{B) >, where Ui{A) and Vp{B) are bases of subspaces 
e{A) and e(-B). That includes creating a 2^^ x 2^^ density matrix p followed by permutations of rows, columns and 
some basic arithmetic operations on the elements of p. Instead of operating on a huge matrix, we pick up only certain 
elements from l^* >, performing basic algebra to build a reduced density matrix directly. Details are in the Appendix. 



III. RESULTS AND DISCUSSIONS 

A. Exact entanglement of a 19-site spin system 

We examine the change of concurrence in Eq. ^ between the center spin and its nearest neighbor as a function of 
A — h/J for both the 7-site and 19-site systems. In Fig[TJ the concurrence of the 7-site system reaches its maximum 
0.15275 when A = 2.61. In the 19-site system, the concurrence reaches 0.0960 when A = 3.95. The maximum value 
of concurrence in the 19-site model, where each site interacts with six neighbors, is roughly 1/3 of the maximum 
concurrence in the one-dimensional transverse Ising model with size N=201 [17], where it has only two neighbors for 
each site. It is the monogamy [30l [3l| that limits the entanglement shared among the number of neighboring sites. 
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However, entanglement between other nearest neighbors are shghtly different than those between the pairs involving 
the center. Figure[2]shows that the less the number of neighbors of a pair the larger the entanglement. The concurrence 
between the 1st and 2nd spins is greater than that between the 1st and 4th in the 7-site system. For the 19-site one, 
the concurrence between the 1st and 2nd spins is greater than that between the 2nd and 5th. The same rule applies 
to the others, therefore €2,5 > C^^ > Cs.io- Although 5 & 6 and 5 & 10 have the same number of neighbors, the 
number of neighbors of neighbors of 5 & 6 is less than that of 5 & 10. Consequently, 6*5,6 is slightly larger than 6*5,10. 

Our numerical calculation shows that the maximum concurrence of next nearest neighbor (say the 1st and 10th 
spins) is less than 10~^. The truly non-local quantum part of the two-point correlations is nonetheless very short- 
ranged Oj]. It shows that the entanglement is short ranged, though global. These results are similar to those obtained 
for Ising one-dimensional spin systems in a transverse magnetic field ^9|]. The range of entanglement, that is the 
maximum distance between two spins at which the concurrence is different from zero, is short. The concurrence 
vanishes unless the two sites are at most next-neighbours. 

B. Entanglement and quantum phase transition 

As we mentioned in Part II-B, all the information needed for quantifying the entanglement of two spins is contained 
in the reduced density matrix obtained from density matrix. In other words, entanglement is coded by the information 
of ground state, while the quantum phase transition is characterized by the change of ground state. In order to 
quantify the change of the many-body wavefunction when the system crosses the critical point, we calculate the 
change of concurrence dC/dX between the center and its nearest neighbor as a function of parameter A for both 7 
and 19 sites systems, as shown in FiglS] It is known that an infinite system, a system is at the thermodynamic limit, 
is supposed to have a singularity at the critical point of quantum phase transition; for a finite system one still has 
to take finite size effect into consideration. However, in Figl3] both systems show strong tendency of being singular 
at Ac — 1.64 and Ac = 3.01 respectively. Renormalization group method for an infinite triangular system predicts 
critical point at Ac — 4.75784; the same method for a square lattice system at Ac — 2.62975 [33|, while finite-size 
scaling has Ac = 3.044 for square lattice [34|]. Our results show that the tendency to be singular is moving towards 
the infinite critical point as the size increases. For one-dimensional system since the calculations can be done for 
large number of spins, finite size scaling calculations for N ranging from 41 to 401 spins indicate that the derivative of 
the concurrence diverges logarithmically with increasing system size [9] . In our study_we can not perform finite size 
scaling analysis since we do not have enough data points to perform data collapse [3^. Optimization methods [3^ 
and Parallel Tracemin code is under development which will allow us to obtain exact results for larger system. 

To understand this model better, a discussion about the degeneracy in the system and an explanation of the 
energy spectrum is necessary. It is known that the ground state degeneracy of the Heisenberg spin model depends 
on whether the total number of spins is even (singlet) or odd (doublet). For the Ising model with transverse field 
on an infinite ID chain, the ground state in the ferromagnetic (FM) phase is doubly degenerate and is gaped from 
the excitation spectrum by 2J(1 — h/J) [s^]- (Note that, however, this degeneracy is never achieved unless one goes 
to the thermaldynamic limit, regardless number of spins being even or odd.) In our model, the Ising coupling is 
ferromagnetic, as opposed to the spin liquid case with antiferromagnetic coupling, the system is expected to break the 
Z2 symmetry and develop the (Ising) FM order under small transverse field. Further, due to the construction of the 
lattice, it is impossible to have a system that has even number of sites while conserving all the lattice group symmetries. 
So we expect that the same doublet degeneracy remains in 2D as the system goes to the thermaldynamic limit. The 
energy gaps from our numerical results of finite systems are less than 10~^ (Fig. 4), which are well consistent with 
the expectation. The strict doublets in finite systems only happen at h/ J = exactly, when entanglements naturally 
are zero, not entangled at all; no matter which one of the doublet ground state is chosen, it gives the same value of 
entanglement. Otherwise even very small h help distinguish the ground state. Technically we don't have to worry 
about that a different superposition of the ground states gives different values of entanglement. 

The energy separation between the ground state and the first excited state in terms of A clarifies the spectrum of 
the system. Figure S] presents the doubly degeneracy of the ground state and they separate around A = 1.5 for the 
7-site and A ~ 2.5 for the 19-site system. While in the dC/dX vs A graph, both systems show strong tendency of 
being singular at Ac = 1.64 and Ac — 3.01 respectively. Both "separation position" and "singular position" are used 
as an indicator of "critical point" . And we believe using the finite-size scaling both will give the same "critical point" 
of the infinite system. But it seems dC /dX vs A is a better indicator because for the same size system it points out 
a value closer to the expected critical point. This property benefits the finite size scaling method, since less/smaller 
systems may be needed. 
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C. Introducing impurities to tune the entanglement 

We introduce one impurity in the center of the 19-site spin system. The impurity only interacts with nearest 
neighbors in strength J' — {1 + a)J . When the strength increases, the concurrence of any two spins decreases. Then 
we move the impurity to site 5. The concurrence shows the same trend of decreasing, as the impurity strengths go 
up. Details are shown in FiglS] for the concurrence of different pairs with various strength of impurity in the center 
of the 7-site system, while FiglHlshow the results for 19-site system. Fig[7] shows the results with various strength of 
impurity at site 1 of a 7-site system and Figl8]for various strength of impurity at site 5 of a 19-site system. 

The maximum of entanglement is shifted with the increasing of the parameter a. The shift is the result of the 
competition between the spin-spin interaction J and the external transverse magnetic field h. Consider the ideal 
situation of pure infinite system. Without the coupling interaction, all the spins will point along the direction of 
transverse field. While with the absence of transverse field, the ground state is supposed to be 2-fold degenerated, 
either along the positive x direction or the negative x direction. Every spin has six neighbors, so averagely is affected 
by three J, and one h. When the two forces are well-matched in strength, the phase transition occurs. Fig. 3 indicates 
the 19-site system has strong tendency of singularity at A = 3.01, which is consistent and very close with the above 
statement. When the system is finite, the boundary effect (less than six neighbors) will affect the position of the 
maximum of entanglement a little bit as Fig. 2 shows for different pairs. After we introduce the impurity, the balance 
of 3 : 1 is destroyed, so the the maximum of entanglement is shifted quite a bit for different strength of ,]' = {! + a) J. 

The value of the maximum changes with a is because of the monogamy that limits the entanglement shared among 
neighbors. For example, in Fig. 4, the stronger the interactions between 1 & 4 and 2 & 4, i.e. the larger a, the less 1 
& 2 entangle. So the value of maximum goes down for larger a. 

Figure [9] gives a good overview of the change of concurrence for the 7-site system. The large yellow dot stands for 
the impurity and silver dots denote regular spins. Lines connecting two sites represent the entanglement. If the line 
is green, it means the entanglement between two sites increases as the impurity gets "stronger" , and the yellow line 
indicates that the entanglement decreases when the impurity increases. We can explain these phenomena of 7-site 
system as follows. When the impurity interacts more with the neighbor, the pair also entangles more. Since some 
spins are more involved with the impurity, they themselves entangle less. The only exceptions are the next nearest 
neighbors. Thus, entanglement close to the impurity tends to get bigger when J' is greater than J. However, the 
behavior of entanglement between site 5 and 10 in the 19-site system surprisingly goes down as the strength of the 
impurity coupling increases. It is not clear why the behavior is different for the one in the 7-site system and whether 
increasing the system size has any effect. We are planing to increase the size of the system to include the next layer, 
which will bring the system to 37-site, in order to analyze this phenomena. 

All the results above are obtained through sequential computing. In the future to increase the object size under 
consideration, we plan to take advantage of parallel computing. We already have a parallel Tracemin algorithm and 
we are developing a parallel code for computation of the partial trace. This will be useful as we expand our 2D 
systems to larger number of spins in order to perform finite-size scaling for quantum critical parameters. 

In summary, the Tracemin algorithm allowed us to carry out an exact calculation of entanglement in a 19-site 
two dimensional transverse Ising model. We demonstrated for such a class of two-dimensional magnetic systems, 
that entanglement can be controlled and tuned by varying the parameter A in the Hamiltonian and by introducing 
impurities into the systems. 
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By studying the patterns of J2<i j> ^ ® • • • o"; 
rules. 



Appendix A: Applications of trace minimization algorithm 
1. General forms of matrix representation of the Hamiltonian 

• • • cr J (81 • • • / and / (g) • • • erf (g) • • • /, one founds the fohowing 



a. (jf for N spins 

The matrix is 2^ by 2^; it has only 2^ diagonal elements. Elements follow the rules shown in Fig. [TU] 
If one stores these numbers in a vector, and initializes v = (iV), then the new v is the concatenation of the original 
V and the original v with 2 subtracted from each of its elements. We repeat this N times, i.e., 



V 

v~2 



(Al) 




(A2) 
(A3) 

(A4) 



E< 



■ / for N spins 



Since 



exclude each other, for matrix / ( 



■ /, every row/column contains 



1 y 1 ^ 

only one "1", then the matrix owns 2^ "l"s and only "1" in it. If we know the position of "l"s, it turns out that 
we can set a 2-'^ by 1 array "col" to store the column position of "l"s corresponding to the 1st 2^th rows. In 
fact, the non-zero elements can be located by the properties stated below. For clarity, let us number N spins in the 
reverse order as: N-1, N-2, . . . , 0, instead of 1, 2, . . . , N. The string of non-zero elements starts from the first row at: 
1 + 2' + 2-'; with string length 2-'; and number of such strings 2^^^^^. For example. Fig. [TT] shows these rules for a 
scheme of / (g) cr? (g tr? €5 J /. 



Again, because of the exclusion, the positions of non-zero element "1" of /g) • • • 



• crj (g) • • • / are different from 



those of / ( 



■I- SoE<,,>^ 



/ is a 2 by 2 matrix with only 1 and 0. 



After storing array "col", we repeat the algorithm for all the nearest pairs < i,j >, and concatenate "coF's to 



position matrix "c" of ^ 



<i,j> 



■ I. In the next section we apply these rules to our problem. 



2. Specialized matrix multiplication 

Using the diagonal elements array "v" of J^i '^t position matrix of non-zero elements "c" for X]<i j> ^® ■ ■ ■ erf ® 
• • • crjig - • • /, we can generate matrix H, representing the Hamiltonian. However, we only need to compute the result of 
the matrix- vector multiplication H*Y in order to run Tracemin, which is the advantage of Tracemin, and consequently 
do not need to explicitly obtain H. Since matrix-vector multiplication is repeated many times throughout iterations, 
we propose an efficient implementation to speedup its computation specifically for Hamiltonian of Ising model (and 
XY by adding one term). 
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For simplicity, first let Y in H*Y be a vector and J = h = 1 (in general Y is a tall matrix and J 7^ /i 7^ 1). Then 



<i,3> 



(I I I \ 

1 1 1 



Y{2) 



(v{l) 



v{2) 



/ r(c(l,l))+y(c(l,2)) + ... + r(c(l,#o/pairs) \ 
y(c(fc, 1)) + y(c(fc, 2)) + . . . + y(c(fc, jfof pairs) 



\ ( YW \ 

Y{2) 



v{2^)) \Y{2^)j 

( v{\)^Y(\) \ 
u(2) *r(2) 



Vy(c(2^,l))+y(c(2^,2)) + ... + r(c(2^,#o/pairs) / V t;(2^ ) * r(2^) / 



(A5) 



where p# stands for the number of pairs. 

When Y is a matrix, we can treat Y (2^ by p) column by column for X]<i j> '8) • • • crf <8) • • • crj (8 • • • /■ Also, we 
can accelerate the computation by treating every row of Y as a vector and adding these vectors at once. Fig. [12] 
visualized the process. 

Notice that the result of the multiplication of the xth row of X)<i j> '''f '''J (delineated by the red box above) and 
Y, is equivalent to the sum of rows of Y, whose row numbers are the column indecis of non-zero elements' of the xth 
row. Such that we transform a matrix operation to straight forward summation & multiplication of numbers. 



Appendix B: Partial Trace 



All the information needed for quantifying the entanglement of two spins i & j is contained in the reduced density 
matrix /r?(j,j), which can be obtained from global density matrix p = |^)('(/'|, where j?/') is the ground state of the 
system, via partial trace. Now let us show how we can obtain the reduced density matrix from the ground state 
calculated by Tracemin. 



1. Density operator in the pure case and partial trace 



Consider a system whose state vector at the instant t is 

n n 

The density operator p(t) is defined as 

p(t) = \i,(t))m)\- (B2) 

It enables us to obtain all the physical predictions of an observable A(t) by 

(A(t)) = Tr{p(t)A}. (B3) 

Let us consider two different system (1) and (2) and the global system (l)-|-(2), whose state space is the tensor 
product: e = e{l) £(2). Let |wn(l)) be a basis of e(l) and |up(2)), a basis of £(2), the kets |u„(l)) |wp(2)) from a 
basis of £. 

The density operator p of the global system is an operator which acts in s. We construct from p an operator p(l) 
[or p{2)] acting only in £(1) [or £(2)] which will enable us to make all the physical prediction about measurements 
bearing only on system(l) or system(2). This operation will be called a partial trace with respect to (2) [or (1)]. 
Matrix elements of the operator p{l) are 

(u„(l)|p(l)|u„-(l)) = 5]((u„(l)|(z;,(2)|)p(|^„,(l))l«p(2))). (B4) 
p 
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Now let A(l) be an observable acting in e{l) and A{1) = A{1) (E) I{2), its extension in e. We obtain, using the 
definition of the trace and closure relation: 

1(1)) = Tr{p{l)A{l)}. (B5) 

As it is designed, the partial trace enables us to calculate all the mean values {A{1)) as if the system(l) were 
isolated and had p{l) for a density operator [s^ . 



2. Properties of the reduced density matrix 

As we calculate the entanglement of formation, we trace out all spins but two (bipartite). Their reduced density 
matrix is, therefore, four by four. Reality and parity conservation of H together with translational invariance already 
fix the structure of p to be symmetric with pn, P22, P33, P44, P14, P23 as the only non-zero entries. It follows from the 
symmetry properties of the Hamiltonian, the p must be real and symmetrical, plus the global phase flip symmetry of 
Hamiltonian, which implies that [a^a^,pij] — 0, so 



1 



/ Pii P12 Pl3 Pl4 

P21 P22 P23 P24 

P3I P32 P33 P3i 

\ P41 Pi2 P43 P44 



/ Pll P12 Pl3 Pli 

P2I P22 P23 P24 

P3l P32 P33 P34 

P41 P42 P43 P44 




= 0, 



Pll 


P12 


Pl3 


Pl4 


-P21 


-P22 


-P23 


-P24 


-P31 


-/O32 


~P33 


~P34 


P41 


P42 


P43 


P44 


( 


P12 


Pl3 


' \ 


-P21 





- 


P24 


-P31 





- 


P34 , 


K 


P42 


P43 


/ 




-Pl3 Pl4 

-p23 P24 

-P33 P34 

-P43 P44 



P12 = Pl3 = P2I = P24 = P31 = P34 = 0. 



Because pij is symmetric. 



therefore, 



Pl4 — P41, P23 — P32-, 



= 0, 



(B6) 

(B7) 

(B8) 
(B9) 

(BIO) 



P^0 = 



Pll pi4 \ 

P22 P23 

P32 P33 

Pil P44 / 



(Bll) 



3. Building the reduced density matrix 

The available code of partial trace involves permutating rows/columns of the density matrix p. On one hand, we 
need to avoid generating the huge matrix p (2^^ by 2^^). On the other hand, even if we have p, permutations are 
too costly to be computed. Fortunately, we are able to convert "generate a global density matrix, then partial trace" 
into "get six elements then build a reduced density matrix". These six elements are closely related to the ground 
state which we have already obtained after the application of Tracemin. In fact the structure of the system (the 
Hamiltonian) guarantees that our ground state is a real vector and implicit of time. That makes the calculation 
neater, with no worry about the complex conjugate and time evolution. 

Here is how to retreat the six elements. For example, we intent on tracing three spins: 1st, 3rd & 5th, out of total 
five, provided the ground state with above properties. We start from the definition of partial trace (IB4p 

(7.„(2)|(t;p(4)|p(2,4)|u„,(2))|v(4)) 
= («a(l)|(^i„(2)|(6/3(3)|(«p(4)|(c^(5)|V')(V'|aa(l))|u„.(2))|6;3(3))|v(4))|c^(5)), (B12) 

a, ^,7 
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and carry it out for a specific matrix element 



Pii(2,4) 
= (00|p(2,4)|00) 

= (a^(l) 6/j(3) c^(5)|V.)(V;|a„(l) 6^3(3) c^(5)) 

a, ^,7 

".■ l^p) is real 

= ^ (Ml) 6^3) Oc^lV)'. (B13) 

a,/3,7 



Since 



IV') = V'm|aa(l)un(2)5^(3)«p(4)c^(5)), (B14) 

(aa(l)u„(2)6/^(3)wp(4)c7(5)|^) = Vm, (B15) 
(00(1) ^0(3) c-y(5)|7/') is the coefBcient in front of base |aQ.(l) bp{3) c-y(5)) at \tp) expansion. When we write 



as a column vector 



/ ^1 \ 100000) 

I jooooi) 

, its elements are the coefficients corresponding to the basis . . Then if we can 

Vv2«y iiiiii) 

locate |aQ(l) 5^(3) c^(5)) as the mth base among all basis, we know that {aa{l) 6/3(3) Cj{5)\ip) — ^m- Our 
task is to locate all the basis with the 2nd & 4th spins being at state |0), then pick out corresponding V'mS, square 
and sum them together. 

Before we continue onto the details, let us construct a basis matrix of five spins illustrated in Fig. [T2] and define 
the following. 



Period : we say the pattern like ^ is a period. 

1 

1 

Segment: : or : is a segment. 
1 

Length: the number of elements in a period or a segment. 
L.P.(i): length of a period of the ith spin. 
L.S.(i): length of a segment of the ith spin. 



The ith spin out of total N spins has: 

ith/N 2nd/h Ath/h 
L.P. 2^-*+i 2^-2+1 ^ iQ 25-4+1 ^ 4 

L.S. 2^-' 2^-2 ^ 8 2^-* ^ 2 

#ofperiods 2^-1 2^-1 = 2 2^-1 = 8 
#of segments 2'' 2^ = 4 2-* = 16 

Using these definitions, we can in general easily locate basis such as | . . . . . . . . .), and within five steps obtain 
element pii(i, j), by applying the following algorithm. 

1. 1st in every period for ith spin is located at "p"= 1, 1 + L.P.{i), 1 + 2L.P.{i), . . . < 2^ 

2. 1st in every period for jth spin, when the ith is 0, is located at "q"= p, p + L.P.{j), 1 + 2L.P.{i), . . . < 
p + L.P.{i) ~ 1, (•.■ I < J .-. L.P.{j) < L.P.ii) kL.S.U) < L.P.ii)). 

3. L.S.(j) decides the length of continued | . . . . . . . . .) basis. 
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4. After locating "q", we naturally have V(9) V'^C?) =^ V'^C?) + V'^C? + 1) + • • • + L.S.{j) - 1), then locate 

the next "q". 

5. When wc add them altogether, it is pii{i, j). 
Similarly, we can locate "01" "10" "11" for 

P22{i,j) = (01|p(i,j)|01), (B16) 

P33(*,j) = (10|p(z,j)|10), (B17) 

P44(«,j) = (ll|p(i,j)|ll). (B18) 

pii{i,j) and P23{i,j) are a bit different. 

Pi4{i,j) = (00|p(j, j)|ll) = (. . . . . . . . . . • . 1 . • . 1 • • •)• (B19) 

In this case we don't have to locate | . . . . . . . . .) &:| . . . 1 . . . 1 . . .) respectively. Letting q be the position of ... ... .. . 
and the corresponding of . . . 1 . . . 1 . . . (i.e. other basis are the same), they are related hy q' = q + L.S.{i) + L.S.{j). 
That enable us obtain pn and pu at the same time: 

Piiihj) = (^20) 

pi4(i,i) = + + (B21) 

So are the P22{i,j) and P23{i,j)- 
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0.18 




FIG. 1: (Color online) Concurrence of center spin and its nearest neighbor as a function of A for both 7-site and 19-site system. 
In the 7-site system, concurrence reaches maximum 0.15275 when A — 2.61. In the 19-site system, concurrence reaches the 
maximum 0.0960 when A = 3.95. 
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FIG. 2: (Color online) The nearest neighbor concurrence as a function of A for different pairs. In 7-site system, there are two 
distinct pairs 1 & 2 and 1 & 4. In 19-site system, they are 1&2, 2&5, 5&6 and 5 & 10. 




FIG. 3: (Color online) The change of concurrence between the center and its nearest neighbor as a function of parameter A for 
7-site and 19-site systems respectively. They both show the strong tendency of singularity at A = 1.54 and A = 3.01. 
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FIG. 4; (Color online) The energy separation between the ground state and the first excited state as a function of A. 
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FIG. 5: (Color online) Concurrences of different pairs with various strength of impurity in the center of the 7-site system. 
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FIG. 7: (Color online) Concurrence of different pairs with various strength of impurity at site 1 of a 7-site system. 
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FIG. 8: (Color online) Concurrence of different pairs with various strength of impurity at site 5 of a 19-site system. 




FIG. 9: (Color online) Overview of the change of concurrence in 7-site system. The large yellow dot stands for the impurity 
and silver dots denote regular spins. Lines connecting two sites represent the entanglement. If the line is green, it means 
the entanglement between two sites increases as the impurity gets "stronger", otherwise the yellow line indicates that the 
entanglement decrease when the impurity increases. 
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FIG. 10: Diagonal elements ol af for N spins. 




FIG. 11; Scheme of matrix / (gjaf^af®/®/ 
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FIG. 12; (Color online) The illustration of H*Y. 
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FIG. 13: (Color online) The illustration of five spins basis. 



